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Abstract: The S— matrix in gravitational high energy scattering is computed from 
the region of large impact parameters b down to the regime where classical gravita- 
tional collapse is expected to occur. By solving the equation of an effective action 
introduced by Amati, Ciafaloni and Veneziano we find that the perturbative ex- 
pansion around the leading eikonal result diverges at a critical value signalling the 
onset of a new regime. We then discuss the main features of our explicitly unitary 
S— matrix down to the Schwarzschild's radius R = 2Ga/s, where it diverges at a 
critical value b ~ 2.25 R of the impact parameter. The nature of the singularity is 
studied with particular attention to the scaling behaviour of various observables at 
the transition. The numerical approach is validated by reproducing the known exact 
solution in the axially symmetric case to high accuracy. 
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1. Introduction 

We discuss in this paper the high-energy quantum string-gravity scattering by using 
the string formulation of Amati, Ciafaloni and Veneziano [1, 2]. Here we concentrate 
on the study of the high-energy elastic ^-matrix for two colliding strings in the 
limit in which the string length X s is negligible with respect both to the gravitational 
Schwarzschild radius R and to the impact parameter b of the the two colliding strings 

^^0, ^^0, R = 2Gy/s, \ s = Va i h, (1.1) 
b b 

where \fs is the center of mass energy. In this point-like-limit (1 = 4 we rely on the 
observation that, because of the soft multi-loop strings amplitudes [1], the S'-matrix 
exponentiates, in terms of an eikonal function of order Gs/H, times a function of 
R 2 /b 2 . In this limit the two colliding strings can be approximated by two massless 
particles 

S = e^ Acl , A d = 2ivGsA d . (1.2) 

Here A c \ is the classical value of an effective action given by a functional of gravita- 
tional fields constructed in terms of effective graviton vertices derived by Lipatov [3] 
in the study of the multi-Regge limit in (QCD and) gravity. The fields involved in the 
effective action A are the longitudinal and the transverse components of the metric 
tensor. The two colliding particles are described in terms of the energy-momentum 
tensor along the two light cone directions x ± = x° ± x 3 . 
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The classical region to study is the case when R is small and it approaches a 
critical value R CT [ t at a fixed values of b. Here the critical value of R is obtained 
by two approximations: we neglect rescattering terms and we neglect the infrared 
singular term of graviton interaction. Within these two approximations one is left 
with the critical behaviour 

R 2 ^ AiO, s) ~ (1 - R/R CTit ) 1/2 Ai ~ (1 - R/R^tf' 2 + constant . (1.3) 
This critical behaviour has been derived in a modified model exhibiting axial sym- 

— * 

metry (where the singular 6— dipole disappears) but here it will be reproduced in the 

— * 

actual case of a dipole at a distance b. 

The paper is organized as follows. First we introduce the model (sec. 2) and 
explain the two relevant approximations for the dipole picture with impact parameter 

— * 

b. Then we explain the solution for R < R cr i t and discuss the search of the critical 
index (sec. 3). We then discuss the results and the behaviour of several observables 
near critical region with R < R cri t- Here we discuss also the region R > -Rcrit in 
the axial case, which is not solvable by our iteration technique. We end with some 
conclusions. 



2. The model 

Here we recall the essential points of the effective action introduced in [1, 2] which are 
needed to describe our numerical study and we describe the two basic simplifications 
that are performed to give 



A 



27rGs J d 2 x [a + t + + a_r - \da + ■ <9a_ + \(tvR) 2 (-(<9 2 0) 2 + 20S a+a _) ) 
S f:9 (x) = d 2 f(x) d 2 g(x) - didjfix) d t d j9 (x) . 

(2.1) 



Here a±(x) are the longitudinal graviton components (dotted lines in Fig. 1). Upon 
factorization the light-cone delta-function 5(x T ) respectively, they depend only on 
the transverse variables x = (xi, x?). In the same way t±(x) correspond to the longi- 
tudinal components of the energy-momentum tensor associated to the two colliding 

— * 

particles at impact parameter b 

t ± (x)=5 2 (x T 1 S- (2-2) 

The field 4>(x) corresponds (upon factorizing the light-cone function 9{x + x )) to the 
transverse graviton component (the wavy line in Fig. la). 

The first two terms of A correspond to the eikonal exchange of longitudinal 
gravitons (Coulomb scattering). This part of the action (R = 0) gives rise to the 
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k = ki + k 2 
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Fig. 1: (a) the scattering vertex between the two partons at distance b. The crosses are 
the external particles. Dotted lines are a± and wiggly lines denote <j>. (b) The rescattering 
of partons. 
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with A an infrared cutoff. The corresponding classical action is given by 



2ttGs 



— ln(A6) 

7T 



(2.3) 



(2.4) 



which shows that the infrared cutoff A is a divergent Coulomb phase contribution 
which does not affect observables. It reproduces Einstein's deflection angle sin = 
R/b. 

The other pieces of A correspond to quantum correction proportional to R 2 /b 2 . 
In particular the last term corresponds to the effective vertex [3] in Fig. la. The 
form of A in Eq. (2.1) is based on two simplifications. 



No rescattering: the effective vertex in Fig. lb, the rescattering term, is ne- 
glected. In this case the dependence of the fields on the light-cone directions 
x + and x~ of the two colliding particles can be factored into the two delta func- 
tions Slx^). Dropping this "trivial" dependence, one is left with the fields in 
the effective action which depend only on the two transverse space components 
x. 
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• No infrared-singular transverse graviton component. The Lipatov vertex of 
Fig. la is given by (9 12 = 9^) 

V£[ T (ki,k 2 ,k) ~ sin 2 9i 2 e^ u T (k) -sin#i 2 cos 9 l2 e^(kj) , (2.5) 
k 2 ^ ' ' 

with €tt and €lt the two independent transverse-traceless tensors of the gravi- 
ton (see for example Eq. 2.5 of Ref. [2]). For k — > the first term sin 2 # 12 //c 2 
remains finite while the second one, sin# 12 cos 9\ 2 /k 2 diverges. In our study 
we neglect the 6lt transverse graviton component which brings in infrared di- 
vergences. The possible role of these singular components and the cancellation 
of the infrared singularities in the ^-matrix has been discussed in [2] by using 
coherent states similar to what is done in QED and QCD [4]. 

In the space representation the first term of the vertex (2.5) gives rise to the 
2<t> S a+ta _ term in the action. Indeed 

F(h)G(k 2 ) sin 2 12 S ft9 (g), (2.6) 

— * 

with f(x) and g{x) the Fourier transforms of the two functions F(k\) and 
G{k 2 ). The second term, infrared singular, gives (ei 2 = e 2 i = — 1) 

F(k 1 )G(k 2 ) sin#i 2 cos#i 2 ddj(x) ■ dd jg (x) e {j (2.7) 

which for the free solution (2.3), with / = a+^ cl , g = a^ cl , is proportional 
to x A b. The contribution from this infrared singular vertex is neglected in 
Eq. (2.1). 

The equations of motion (2.1) are readily obtained by using the fact that (f)S a+a _ = 
a± S^ a ^ up to total derivatives. One has 

d 2 a ± (x) = -2 5 2 (x ± \b) + 2(irR) 2 S^ a± (x) 

(d 2 ) 2 <t>(x) = -S a+a _(x) (2 ' 8) 

If a±(x),(p(x) is a solution, then 

a±(x) = a T {b — x) , (f)(x) = <j){b — x) . (2.9) 

is also a solution. If the solution is unique (e.g. for small R) then it is symmetric. 
Perturbative corrections in R 2 /b 2 to Eq. (2.3) are finite. To first order the action is 
given by 

^>- 1 /-,„(A6) + S + ...Y (2.10) 



2nGs tt V 8b 2 

Here one has used (2.1) thus neglecting the transverse e LT graviton component. This 
latter has been computed in [2] and gives an additional contribution R 2 /(8irb 2 ). 
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3. Solution by an iterative algorithm 

In the small R regime we solve Eq.(2.8) by iteration (A = d 2 ) 

a±(x) <- (1 - u)a±(x) - u A" 1 (2 5 2 {x t 6/2) - 2(n R) 2 ^, a± ) (3.1) 
0(f) <- (1 - cu)0(f ) - A- 2 S a+ , a _ (3.2) 

(5*/^ as given in Eq. (2.1)). Here u is a "relaxation parameter" which may help 
convergence; in a linear context, the iteration converges if the r.h.s. has norm less 
than one, using the relaxation parameter can be useful if all eigenvalues have real part 
less than one. In our non-linear context the role of relaxation is not immediately clear. 
It is only used to check that divergence of the iteration is not changed by a trivial 
modification of the iteration mechanism. All results presented later on correspond 
to u — 1. 

In the actual implementation of the algorithm, we substitute the Dirac deltas 
with extended normalized Gaussian source terms 

5 2 (x ± 6/2) -ATexpj ~ (f ^ /2)2 | (3.3) 

since dealing with smooth functions improves the numerical stability; the Gaussian 
width a enters as a smearing parameter in the calculation. 

3.1 Algorithmic details 

The calculation is organized as follows. We choose a finite 2— dimensional lattice 
of size 2L and lattice spacing 2L/N, typically L = 64, A = 256; we introduce the 
dual lattice (fci,/^) and define the basic differential operators in terms of ki, i.e. 
A — > —k\ — k 2 , etc. The r.h.s. of the iteration map is then computed by going to 
Fourier space where necessary and taking local products of fields in coordinate space. 
This is done very efficiently for any N by using FFTW [5], the best implementation 
of the fast Fourier transform to our knowledge. The main problem is to monitor the 
iteration and stop it when it can be decided that we have convergence or divergence. 
To monitor the iteration we choose a specific component of the action, namely the 
part proportional to R 2 (A = 2ttGs A) 

M = ±(nR) 2 J d 2 x (-(9 2 0) 2 + 20S a+a _) ^A^+A 4>a+a _=R 2 ^A. (3.4) 

A second monitoring device is provided by the fact, already used in [2], that on the 
solutions one has 

2„V + «V+a_=0. (3.5) 

When M shows a variation less that a certain prefixed value (say 1CT 6 ) then we 
stop the iteration and proceed to a more stringent convergence check based on the 



- 5 - 



analysis of the linearized equation near the approximate solution 1 . This is done as 
follows: let 

/C(a ± ,0) (3.6) 

denote the iteration map; near convergence let <fi = <fi + 5<p, a± = a± + 5a± where 
0, a± represent the approximate solution obtained at the last iteration. Then we 
have 

ft) - (£; S ft) 

where the matrix Dy(a±,0) is obtained by linearizing the equation and discarding 
quadratic terms in 5a±, 5(f). The spectrum of is then computed. If the spec- 
trum is contained in the unit circle then we conclude for convergence, otherwise for 
divergence 2 . The linear operator T> contains partial differential operators and the 
fields a±, 0. Its largest eigenvalue, in absolute value, is computed 3 by the Arnoldi 
algorithm contained in the mathematical library ARPACK [6]. 

3.2 The search for _R crit 

Running the algorithm, we observe that the iteration converges at small values of 
R, where we effectively obtain the perturbative solution. However, increasing the 
parameter one finds that the algorithm ceases to converge at a value of R = R crit 
proportional to the impact parameter. 

To identify the value of i? cr it where the transition occurs we adopt a "bisection 
method". Namely we start from an interval [R m i n , R m ax] and set R = \{R m in J rRmax)\ 
if the iteration starting from R is divergent we set R ma x = R otherwise R m in = R\ 
the process is replicated until Rmax — Rmin is less than a desired accuracy (typically 

io- 5 ). 

The value thus obtained is however affected by the presence of "technical" param- 
eters whose impact on the calculation must be carefully analyzed. The parameters 
are L,N,u,a, i.e. infrared and ultraviolet cutoffs, the relaxation parameter and the 
width of the sources. On dimensional grounds one has 

R cri t = bF N ( ( T/L 2 ,b/L,tj), (3.8) 

and, in principle, the limit o — > 0, L — > oo, iV — > oo should be taken. It turns out 
that the iteration scheme is not sensitive to the relaxation parameter, hence we drop 

1 For R > R cr it the iteration will diverge and the program will go in overflow unless we apply 
some strategy to early identify a divergent behaviour. A possible method is to monitor the curvature 
of the function M. as a function of the iteration count. 

2 As an option we may consider the eigenvalue with largest real part: if this is larger than one 
we have divergence, otherwise convergence could be achieved by a suitable choice of relaxation 
parameter. 

3 Recent versions of Matlab© provide a user-friendly interface to the package. 
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it from now on. On the other hand the dependence on infrared (L) and ultraviolet 
(er) parameters is quite accentuated and some fitting procedure must be adopted in 
order to get the continuum infinite volume limit. Dependence on iV on the other 
hand is very weak above N ~ 100, as we argue in next paragraph. 

Continuum limit and Finite Size Scaling. To estimate the N dependence we 
performed a series of iterations starting with different parameters but keeping the 
adimensional ratios inside Fjv fixed. One finds that the results reach a flat plateau 
very soon (above N=100 the numbers agree to 4 figures). This is surely due to 
the way differential operators are dealt with: any finite-difference scheme would 
introduce a systematic error of 0(N~ n ) - e.g. the roughest scheme corresponds 
to d 2 — > 2^-(l — cos(aki))/a 2 ; using k 2 itself in Fourier space makes the error 
exponentially small in N. 

Having settled the N dependence, systematic errors coming from o and L should 
be addressed. A long run on a lattice of values for b, a and L manifests a good scaling 
behaviour (i.e. the plot of R 2 /b 2 is acceptably smooth, hence showing compatibility 
with a scaling behaviour). 

4. Results 

The main result regards the dependence of the critical value of R with respect to 
the impact parameter. We ran the program for a three-dimensional grid of values 
a <C fe <C L) and we looked for a best fit to the F N in Eq. (3.8) assuming a simple 
form 

F N = v/1 + ai a/b 2 + a 2 (b/L) 2 + a 3 (b/L)* + .... (4.1) 
This gives (see Fig. 2) 

i? crit w 0.445(1) b . (4.2) 

The evidence from the fit is that the correct slope corresponds to the envelope of the 
data from below, which can be interpreted in the sense that ultraviolet and infrared 
cutoffs both make the system more stable against collapse. 

The axially symmetric case. A comparison with the result [2, 7] obtained in the 
case of axial symmetry, namely R C nt/b = 2 1 / 2 /3 3 / 4 0.6204 will support our analysis. 
In this case one of the sources is uniformly distributed around an annulus of radius 
b, while the other is localized at the center. The difference between this result and 
the value in Eq. (4.2) suggests that the critical value, due to the non-linearity of 
the equations, may depend on the sources. In order to validate the accuracy of our 
numerical code, we applied it to the study of this axially symmetric case. A fit 
conducted with the same systematics as before gave the result 

R crit ~ 0.6216 (4.3) 
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Fig. 2: i? cr ;t vs. b in the dipole model. 

which presents the same pattern as in Fig. 2 and reproduces the exact slope to an 
accuracy of 0.5%; this check makes us confident on the accuracy of our code and 
allows us to estimate the error in Eq. (4.2) to less than 1%. 

4.1 Critical behaviour 

In order to better understand the nature of the transition at R = R C mt we shall now 
present some results about the critical behaviour of certain observables. The main 
fact we derive from our numerical data is that all observables that we examined have 
a scaling behaviour near the transition which can be reproduced very accurately by 
a square root singularity. This fact supports the conclusion that we are in presence 
of a genuine transition and not simply a breakdown of the iteration scheme. The 
argument is as follows: the iteration scheme represents an efficient way to sum up 
the perturbative expansion in the parameter /C = 2(nR) 2 ; as such the iteration's con- 
vergence radius is regulated by the nearest singularity in the complex fC plane. Our 
analysis shows that the divergence of the iteration scheme is caused by a singularity 
on the real line, which must then correspond to a physical singularity. 

Spectral properties. Let's start with the spectrum of the linearized equation 
which is used in monitoring the convergence of the iteration algorithm. From the 
data a scaling property emerges which appears to be rather robust against variations 
of other parameters, namely the dependence of the spectral radius against R. Let Ao 
denote the spectral radius of the linearized equation: the plot of 1 — Ao as a function 
of a/1 — R/Rcvit is reported in the next picture (Fig. 3) and it suggests a relation of 
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Fig. 3: log-log plot of the spectral radius as a function of R in the dipole model. 



with C close to one. 

The plot reports data from different values of b and it hints at the fact that 
there is universality in the relation Ao = Ao(-R/-R C rit), i-e. the dependence on the 
impact parameter is only through the function R cr it(b) and it is quite insensitive to 
the various cutoff parameters. 

A further investigation regards the region R > i? crit in the case of axially sym- 
metric sources. It is known [2, 7] that above the transition the solution becomes 
complex. We computed the spectrum of the linearized equation around the exact so- 
lution. The spectral radius is compatible with the analytic continuation of Eq. (4.4), 
namely the largest eigenvalue gets an imaginary part - see Fig. 5. 

Classical Action, hrr an d ^lt- Near criticality we examined other observables, 
namely global quantities like 

• the action piece M. as in Eq. (3.4); 

• J cPx h 2 LT , which does not enter the action in the present approximation; 

• hj"j- = A(j), integrated over a finite region oc b. 

As expected, we observed the same square root behaviour for M. 



the type 




(4.4) 




(4.5) 



and for other cases - see Fig. 6. By integrating over R, one gets Eq. (1.3). 
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Fig. 4: The scaling of the leading eigen- Fig. 5: The leading eigenvalue above 
value in the axial symmetric case below -Rcrit in the axial symmetric case above 
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Fig. 6: Scaling behaviour of Ai, Htt and /ilt near Local observables. 



5. Conclusions 

We have neglected the string-size effects for R , b ^> A s (with X s the string effect) and 
studied the region without rescattering and without trace infrared effects for R small 
compared to b in (2.1). We arrived at the critical region in (4.5) with the square root 
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behaviour ensured for various decades. This result nicely agrees with that obtained 
in the axially symmetric case [2, 7]. The fact that at the transition an eigenvalue 1 
appears in the spectrum of the linearized equation shows that the solution (if any) 
at R > R crit is unstable against small perturbations, hence features that have been 
considered irrelevant in the derivation of Eq. (2.1) may come back as relevant. 
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